{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mesh generation for tutorial 06, case 7a\n",
    "\n",
    "This file generates the mesh which is used in the following examples:\n",
    "* 7a_stokes_dirichlet_control\n",
    "\n",
    "The test case is from section 5 of\n",
    "```\n",
    "F. Negri, A. Manzoni and G. Rozza. Reduced basis approximation of parametrized optimal flow control problems for the Stokes equations. Computer and Mathematics with Applications, 69(4):319-336, 2015.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import gmsh\n",
    "from mpi4py import MPI\n",
    "from dolfinx.io import XDMFFile\n",
    "from multiphenicsx.mesh.utils import gmsh_to_fenicsx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert MPI.COMM_WORLD.size == 1, \"This mesh generation notebook is supposed to be run in serial\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Geometrical parameters and related quantities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L1 = 0.9\n",
    "L2 = 0.35\n",
    "L3 = 0.55\n",
    "L4 = 0.2\n",
    "H = 1.0\n",
    "r = 0.1\n",
    "lcar = 0.025"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create mesh with gmsh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()\n",
    "gmsh.model.add(\"mesh\")\n",
    "p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, lcar)\n",
    "p1 = gmsh.model.geo.addPoint(L1, 0.0, 0.0, lcar)\n",
    "p2 = gmsh.model.geo.addPoint(L1 + L2, 0.0, 0.0, lcar)\n",
    "p3 = gmsh.model.geo.addPoint(L1 + L2 + L3, 0.0, 0.0, lcar)\n",
    "p4 = gmsh.model.geo.addPoint(L1 + L2 + L3 + L4, 0.0, 0.0, lcar)\n",
    "p5 = gmsh.model.geo.addPoint(L1 + L2 + L3 + L4, H, 0.0, lcar)\n",
    "p6 = gmsh.model.geo.addPoint(L1 + L2 + L3, H, 0.0, lcar)\n",
    "p7 = gmsh.model.geo.addPoint(L1 + L2, H, 0.0, lcar)\n",
    "p8 = gmsh.model.geo.addPoint(L1, H, 0.0, lcar)\n",
    "p9 = gmsh.model.geo.addPoint(0.0, H, 0.0, lcar)\n",
    "p10 = gmsh.model.geo.addPoint(L1, H / 2, 0.0, lcar)\n",
    "p11 = gmsh.model.geo.addPoint(L1, H / 2 + r, 0.0, lcar)\n",
    "p12 = gmsh.model.geo.addPoint(L1, H / 2 - r, 0.0, lcar)\n",
    "p13 = gmsh.model.geo.addPoint(L1 + L2, H / 2 - r, 0.0, lcar)\n",
    "p14 = gmsh.model.geo.addPoint(L1 + L2 + L3, H / 2 - 3 * r, 0.0, lcar)\n",
    "p15 = gmsh.model.geo.addPoint(L1 + L2 + L3, H / 2 + 3 * r, 0.0, lcar)\n",
    "p16 = gmsh.model.geo.addPoint(L1 + L2, H / 2 + r, 0.0, lcar)\n",
    "l0 = gmsh.model.geo.addLine(p0, p1)\n",
    "l1 = gmsh.model.geo.addLine(p1, p2)\n",
    "l2 = gmsh.model.geo.addLine(p2, p3)\n",
    "l3 = gmsh.model.geo.addLine(p3, p4)\n",
    "l4 = gmsh.model.geo.addLine(p4, p5)\n",
    "l5 = gmsh.model.geo.addLine(p5, p6)\n",
    "l6 = gmsh.model.geo.addLine(p6, p7)\n",
    "l7 = gmsh.model.geo.addLine(p7, p8)\n",
    "l8 = gmsh.model.geo.addLine(p8, p9)\n",
    "l9 = gmsh.model.geo.addLine(p9, p0)\n",
    "l10 = gmsh.model.geo.addLine(p12, p13)\n",
    "l11 = gmsh.model.geo.addLine(p13, p14)\n",
    "l12 = gmsh.model.geo.addLine(p14, p15)\n",
    "l13 = gmsh.model.geo.addLine(p15, p16)\n",
    "l14 = gmsh.model.geo.addLine(p16, p11)\n",
    "l15 = gmsh.model.geo.addLine(p13, p16)\n",
    "l16 = gmsh.model.geo.addLine(p1, p12)\n",
    "l17 = gmsh.model.geo.addLine(p11, p8)\n",
    "l18 = gmsh.model.geo.addLine(p2, p13)\n",
    "l19 = gmsh.model.geo.addLine(p16, p7)\n",
    "l20 = gmsh.model.geo.addLine(p3, p14)\n",
    "l21 = gmsh.model.geo.addLine(p15, p6)\n",
    "c0 = gmsh.model.geo.addCircleArc(p11, p10, p12)\n",
    "line_loop_subdomain1 = gmsh.model.geo.addCurveLoop([l0, l16, -c0, l17, l8, l9])\n",
    "line_loop_subdomain2a = gmsh.model.geo.addCurveLoop([l1, l18, -l10, -l16])\n",
    "line_loop_subdomain2b = gmsh.model.geo.addCurveLoop([l7, -l17, -l14, l19])\n",
    "line_loop_subdomain3a = gmsh.model.geo.addCurveLoop([l2, l20, -l11, -l18])\n",
    "line_loop_subdomain3b = gmsh.model.geo.addCurveLoop([l6, -l19, -l13, l21])\n",
    "line_loop_subdomain3c = gmsh.model.geo.addCurveLoop([l3, l4, l5, -l21, -l12, -l20])\n",
    "line_loop_subdomain4 = gmsh.model.geo.addCurveLoop([l11, l12, l13, -l15])\n",
    "subdomain1 = gmsh.model.geo.addPlaneSurface([line_loop_subdomain1])\n",
    "subdomain2a = gmsh.model.geo.addPlaneSurface([line_loop_subdomain2a])\n",
    "subdomain2b = gmsh.model.geo.addPlaneSurface([line_loop_subdomain2b])\n",
    "subdomain3a = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3a])\n",
    "subdomain3b = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3b])\n",
    "subdomain3c = gmsh.model.geo.addPlaneSurface([line_loop_subdomain3c])\n",
    "subdomain4 = gmsh.model.geo.addPlaneSurface([line_loop_subdomain4])\n",
    "gmsh.model.geo.synchronize()\n",
    "gmsh.model.addPhysicalGroup(1, [l9], 1)\n",
    "gmsh.model.addPhysicalGroup(1, [l0, l1, l2, l3, l5, l6, l7, l8], 2)\n",
    "gmsh.model.addPhysicalGroup(1, [l4], 3)\n",
    "gmsh.model.addPhysicalGroup(1, [l10, l14], 4)\n",
    "gmsh.model.addPhysicalGroup(1, [c0, l15], 5)\n",
    "gmsh.model.addPhysicalGroup(2, [subdomain1], 1)\n",
    "gmsh.model.addPhysicalGroup(2, [subdomain2a, subdomain2b], 2)\n",
    "gmsh.model.addPhysicalGroup(2, [subdomain3a, subdomain3b, subdomain3c], 3)\n",
    "gmsh.model.addPhysicalGroup(2, [subdomain4], 4)\n",
    "gmsh.model.mesh.generate(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert to a dolfinx mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, subdomains, boundaries = gmsh_to_fenicsx(gmsh.model, gdim=2)\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save mesh, subdomains and boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with XDMFFile(MPI.COMM_WORLD, \"vorticity_reduction.xdmf\", \"w\") as file:\n",
    "    file.write_mesh(mesh)\n",
    "    file.write_meshtags(subdomains)\n",
    "    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)\n",
    "    file.write_meshtags(boundaries)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
